Exon

Grange: class. includes genomic locations & associations * seqnames: could be chromosomes, contigs, … * star (*) means unknown strand

# + = plus strand 
exon = GRanges(c('chr1:20-30:+', 'chr1:40-50:+', "chr1:45-55")) 
exon

## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1     20-30      +
##   [2]     chr1     40-50      +
##   [3]     chr1     45-55      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

exon = GRanges(c('chr1:20-30:+', 'chr1:40-50:+', "chr1:45-55:+")) 

exon is a vector (like DNAstring)

length(exon)

## [1] 3

start coordinate

start(exon)

## [1] 20 40 45

end(exon)

## [1] 30 50 55

Implementation of bioconductor assumes that intervals are closed intervals (in the operations) * e.g. 40 to 50: 50 - 40 + 1 = 11

width(exon)

## [1] 11 11 11

Collapse all exons regions

reduce(exon)

## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1     20-30      +
##   [2]     chr1     40-55      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Disjoint range

disjoin(exon)

## GRanges object with 4 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1     20-30      +
##   [2]     chr1     40-44      +
##   [3]     chr1     45-50      +
##   [4]     chr1     51-55      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

promoters(exon, upstream=2000, downstream=3000)

## GRanges object with 3 ranges and 0 metadata columns:
##       seqnames     ranges strand
##          <Rle>  <IRanges>  <Rle>
##   [1]     chr1 -1980-3019      +
##   [2]     chr1 -1960-3039      +
##   [3]     chr1 -1955-3044      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

Project reads onto the genome * Interpretation: 19 zeros, 11 ones, 9 zeros, 5 ones, 6 twos, 5 ones * (see representation in next chunk for mvisual representation)

coverage(exon)

## RleList of length 1
## $chr1
## integer-Rle of length 55 with 6 runs
##   Lengths: 19 11  9  5  6  5
##   Values :  0  1  0  1  2  1

coverage(exon) %>% unlist() %>% as.integer()

##  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0
## [36] 0 0 0 0 1 1 1 1 1 2 2 2 2 2 2 1 1 1 1 1

Coerce it as granges

cvg <- coverage(exon)
cvg <- cvg %>% as('GRanges')
cvg

## GRanges object with 6 ranges and 1 metadata column:
##       seqnames    ranges strand |     score
##          <Rle> <IRanges>  <Rle> | <integer>
##   [1]     chr1      1-19      * |         0
##   [2]     chr1     20-30      * |         1
##   [3]     chr1     31-39      * |         0
##   [4]     chr1     40-44      * |         1
##   [5]     chr1     45-50      * |         2
##   [6]     chr1     51-55      * |         1
##   -------
##   seqinfo: 1 sequence from an unspecified genome

granges(cvg)

## GRanges object with 6 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1      1-19      *
##   [2]     chr1     20-30      *
##   [3]     chr1     31-39      *
##   [4]     chr1     40-44      *
##   [5]     chr1     45-50      *
##   [6]     chr1     51-55      *
##   -------
##   seqinfo: 1 sequence from an unspecified genome

You can add elements to this dataframe

cvg$score <- c(1,0,0,4,3,5)
cvg

## GRanges object with 6 ranges and 1 metadata column:
##       seqnames    ranges strand |     score
##          <Rle> <IRanges>  <Rle> | <numeric>
##   [1]     chr1      1-19      * |         1
##   [2]     chr1     20-30      * |         0
##   [3]     chr1     31-39      * |         0
##   [4]     chr1     40-44      * |         4
##   [5]     chr1     45-50      * |         3
##   [6]     chr1     51-55      * |         5
##   -------
##   seqinfo: 1 sequence from an unspecified genome

mcols(cvg)

## DataFrame with 6 rows and 1 column
##       score
##   <numeric>
## 1         1
## 2         0
## 3         0
## 4         4
## 5         3
## 6         5

df <- DataFrame(
  i   = 1:3,
  dna = DNAStringSet(c('AAA','CCC','GGG')),
  gr  = exon
)
df

## DataFrame with 3 rows and 3 columns
##           i            dna           gr
##   <integer> <DNAStringSet>    <GRanges>
## 1         1            AAA chr1:20-30:+
## 2         2            CCC chr1:40-50:+
## 3         3            GGG chr1:45-55:+

SNP

snp <- GRanges(c('chr1:23435','chr2:54323'))

Shift SNP locations (e.g. if another reference was used):

shift(snp, 1)

## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1     23436      *
##   [2]     chr2     54324      *
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

Extend SNP ranges:

flank(snp, 10)

## GRanges object with 2 ranges and 0 metadata columns:
##       seqnames      ranges strand
##          <Rle>   <IRanges>  <Rle>
##   [1]     chr1 23425-23434      *
##   [2]     chr2 54313-54322      *
##   -------
##   seqinfo: 2 sequences from an unspecified genome; no seqlengths

How many of the exons contain SNPs inside, what is the overlap?

snp <- GRanges(c('chr1:15:+','chr1:25:+','chr1:35:+'))
countOverlaps(snp, exon)

## [1] 0 1 0

Viceversa:

countOverlaps(exon, snp)

## [1] 1 0 0

findOverlaps(snp, exon)

## Hits object with 1 hit and 0 metadata columns:
##       queryHits subjectHits
##       <integer>   <integer>
##   [1]         2           1
##   -------
##   queryLength: 3 / subjectLength: 3

snp[snp %over% exon]

## GRanges object with 1 range and 0 metadata columns:
##       seqnames    ranges strand
##          <Rle> <IRanges>  <Rle>
##   [1]     chr1        25      +
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

sessionInfo()

## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.3
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] GenomicRanges_1.36.0 GenomeInfoDb_1.20.0  ggbeeswarm_0.6.0    
##  [4] forcats_0.4.0        stringr_1.4.0        dplyr_0.8.2         
##  [7] purrr_0.3.2          readr_1.3.1          tidyr_0.8.3         
## [10] tibble_2.1.3         ggplot2_3.2.0        tidyverse_1.2.1     
## [13] Biostrings_2.52.0    XVector_0.24.0       IRanges_2.18.1      
## [16] S4Vectors_0.22.0     BiocGenerics_0.30.0 
## 
## loaded via a namespace (and not attached):
##  [1] beeswarm_0.2.3         tidyselect_0.2.5       xfun_0.8              
##  [4] haven_2.1.0            lattice_0.20-38        colorspace_1.4-1      
##  [7] generics_0.0.2         htmltools_0.3.6        yaml_2.2.0            
## [10] rlang_0.4.0            pillar_1.4.2           glue_1.3.1            
## [13] withr_2.1.2            modelr_0.1.4           readxl_1.3.1          
## [16] GenomeInfoDbData_1.2.1 zlibbioc_1.30.0        munsell_0.5.0         
## [19] gtable_0.3.0           cellranger_1.1.0       rvest_0.3.4           
## [22] evaluate_0.14          knitr_1.23             vipor_0.4.5           
## [25] broom_0.5.2            Rcpp_1.0.1             scales_1.0.0          
## [28] backports_1.1.4        jsonlite_1.6           hms_0.4.2             
## [31] digest_0.6.19          stringi_1.4.3          grid_3.6.0            
## [34] bitops_1.0-6           cli_1.1.0              tools_3.6.0           
## [37] magrittr_1.5           RCurl_1.95-4.12        lazyeval_0.2.2        
## [40] crayon_1.3.4           pkgconfig_2.0.2        xml2_1.2.0            
## [43] lubridate_1.7.4        assertthat_0.2.1       rmarkdown_1.13        
## [46] httr_1.4.0             rstudioapi_0.10        R6_2.4.0              
## [49] nlme_3.1-140           compiler_3.6.0

Beeswarm plots

jtrplot <- p +
    geom_jitter(aes(colour = gene)) +
    theme(legend.position = "none")
jtrplot

dotplot <- p + geom_dotplot(binaxis = "y", binwidth = 1/6,
                            stackdir = "center", stackratio = 0.75,
                            aes(color = gene)) +
    theme(legend.position = "none")
dotplot

Beeplot: to avoid overlapping points

beeplot <- p + geom_beeswarm(aes(color = gene)) + 
    theme(legend.position = "none")
beeplot

# jtrplot + dotplot + beeplot

dfx <- as.data.frame(Biobase::exprs(hiiragi2013))

scp <- ggplot(dfx, aes(x=`59 E4.5 (PE)`, y=`92 E4.5 (FGF4-KO)`)) + geom_point()
scp

scp <- ggplot(dfx, aes(x=`59 E4.5 (PE)`, y=`92 E4.5 (FGF4-KO)`)) + geom_point(alpha=0.1)
scp

scp <- ggplot(dfx, aes(x=`59 E4.5 (PE)`, y=`92 E4.5 (FGF4-KO)`)) + geom_density2d(h=0.5, bins=60)
scp

scp <- ggplot(dfx, aes(x=`59 E4.5 (PE)`, y=`92 E4.5 (FGF4-KO)`)) + geom_density2d()
scp

scp <- ggplot(dfx, aes(x=`59 E4.5 (PE)`, y=`92 E4.5 (FGF4-KO)`)) + geom_hex()
scp

Facet grid

Plotly

Useful to check annotations of data points

library("plotly")

scp <- ggplot(dfx[1:100, ],
              aes(x= `59 E4.5 (PE)`, y = `92 E4.5 (FGF4-KO)`))

scp2 <- scp + geom_point()
ggplotly(scp2)

Exercises

Overall sample counts quantification * Representation of technical variation between samples

colSums(assay(airway))

## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 
##   20637971   18809481   25348649   15163415   24448408   30818215 
## SRR1039520 SRR1039521 
##   19126151   21164133

• Average number of reads per gene (not necessarily gene expression, because read counts depend on gene length, etc.)

genemeans <- rowMeans(log(assay(airway) + 1, base = exp(1)))
ggplot(data=data.frame(genemeans),
       aes(x=genemeans)) +
  geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

org packages

Look for packages with org in their name

• org = organism-based
• hs = homo sapiens
• eg = entrez gene id are central symbols

BiocManager::available("^org\\.")

##  [1] "org.Ag.eg.db"      "org.At.tair.db"    "org.Bt.eg.db"     
##  [4] "org.Ce.eg.db"      "org.Cf.eg.db"      "org.Dm.eg.db"     
##  [7] "org.Dr.eg.db"      "org.EcK12.eg.db"   "org.EcSakai.eg.db"
## [10] "org.Gg.eg.db"      "org.Hs.eg.db"      "org.Mm.eg.db"     
## [13] "org.Mmu.eg.db"     "org.Pf.plasmo.db"  "org.Pt.eg.db"     
## [16] "org.Rn.eg.db"      "org.Sc.sgd.db"     "org.Ss.eg.db"     
## [19] "org.Xl.eg.db"

# BiocManager::install("org.Hs.eg.db")
library(org.Hs.eg.db)
library(AnnotationHub)
org <- org.Hs.eg.db

keytypes availabel for org?

keytypes(org)

##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT" 
##  [5] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"    
##  [9] "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"       
## [13] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"    
## [17] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
## [25] "UNIGENE"      "UNIPROT"

Things that we can map to

columns(org)

##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT" 
##  [5] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"    
##  [9] "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"       
## [13] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"    
## [17] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"      
## [25] "UNIGENE"      "UNIPROT"

Mappings can be done via: * 1 to 1 mapping * 1 to many

ensid <- sample(rownames(airwaynonull), 10)

Mappings via mapIds

mapIds(org, ensid, "SYMBOL", "ENSEMBL") %>% tibble::enframe("Ensemble","Symbol")

## 'select()' returned 1:1 mapping between keys and columns

## # A tibble: 10 x 2
##    Ensemble        Symbol   
##    <chr>           <chr>    
##  1 ENSG00000267086 <NA>     
##  2 ENSG00000184260 HIST2H2AC
##  3 ENSG00000189280 GJB5     
##  4 ENSG00000253121 <NA>     
##  5 ENSG00000237857 <NA>     
##  6 ENSG00000154451 GBP5     
##  7 ENSG00000264201 MIR4701  
##  8 ENSG00000259098 <NA>     
##  9 ENSG00000144369 FAM171B  
## 10 ENSG00000180658 OR2A4

mapIds(org, ensid, "GENENAME", "ENSEMBL") %>% tibble::enframe("GeneName","ENSEMBL")

## 'select()' returned 1:1 mapping between keys and columns

## # A tibble: 10 x 2
##    GeneName        ENSEMBL                                         
##    <chr>           <chr>                                           
##  1 ENSG00000267086 <NA>                                            
##  2 ENSG00000184260 histone cluster 2 H2A family member c           
##  3 ENSG00000189280 gap junction protein beta 5                     
##  4 ENSG00000253121 <NA>                                            
##  5 ENSG00000237857 <NA>                                            
##  6 ENSG00000154451 guanylate binding protein 5                     
##  7 ENSG00000264201 microRNA 4701                                   
##  8 ENSG00000259098 <NA>                                            
##  9 ENSG00000144369 family with sequence similarity 171 member B    
## 10 ENSG00000180658 olfactory receptor family 2 subfamily A member 4

Multiple mappings:

AnnotationDbi::select(org, ensid, c("ENTREZID", "SYMBOL"), "ENSEMBL")

## 'select()' returned 1:1 mapping between keys and columns

##            ENSEMBL  ENTREZID    SYMBOL
## 1  ENSG00000267086      <NA>      <NA>
## 2  ENSG00000184260      8338 HIST2H2AC
## 3  ENSG00000189280      2709      GJB5
## 4  ENSG00000253121      <NA>      <NA>
## 5  ENSG00000237857      <NA>      <NA>
## 6  ENSG00000154451    115362      GBP5
## 7  ENSG00000264201 100616262   MIR4701
## 8  ENSG00000259098      <NA>      <NA>
## 9  ENSG00000144369    165215   FAM171B
## 10 ENSG00000180658     79541     OR2A4

AnnotationDbi::select(org, ensid, c("GO"), "ENSEMBL")

## 'select()' returned 1:many mapping between keys and columns

##            ENSEMBL         GO EVIDENCE ONTOLOGY
## 1  ENSG00000267086       <NA>     <NA>     <NA>
## 2  ENSG00000184260 GO:0000786      IEA       CC
## 3  ENSG00000184260 GO:0003674       ND       MF
## 4  ENSG00000184260 GO:0005634      HDA       CC
## 5  ENSG00000184260 GO:0005634      IDA       CC
## 6  ENSG00000184260 GO:0008150       ND       BP
## 7  ENSG00000184260 GO:0046982      IEA       MF
## 8  ENSG00000184260 GO:0070062      HDA       CC
## 9  ENSG00000189280 GO:0005922      IEA       CC
## 10 ENSG00000189280 GO:0007154      IEA       BP
## 11 ENSG00000189280 GO:0008544      TAS       BP
## 12 ENSG00000189280 GO:0016021      IEA       CC
## 13 ENSG00000189280 GO:0060707      IEA       BP
## 14 ENSG00000189280 GO:0060708      IEA       BP
## 15 ENSG00000189280 GO:0060713      IEA       BP
## 16 ENSG00000189280 GO:1905867      IEA       BP
## 17 ENSG00000253121       <NA>     <NA>     <NA>
## 18 ENSG00000237857       <NA>     <NA>     <NA>
## 19 ENSG00000154451 GO:0000139      IEA       CC
## 20 ENSG00000154451 GO:0003924      IEA       MF
## 21 ENSG00000154451 GO:0005515      IPI       MF
## 22 ENSG00000154451 GO:0005525      IEA       MF
## 23 ENSG00000154451 GO:0005737      IDA       CC
## 24 ENSG00000154451 GO:0005794      IDA       CC
## 25 ENSG00000154451 GO:0006954      IEA       BP
## 26 ENSG00000154451 GO:0009617      IEA       BP
## 27 ENSG00000154451 GO:0016020      HDA       CC
## 28 ENSG00000154451 GO:0031410      IEA       CC
## 29 ENSG00000154451 GO:0034067      IDA       BP
## 30 ENSG00000154451 GO:0042802      IPI       MF
## 31 ENSG00000154451 GO:0042803      IDA       MF
## 32 ENSG00000154451 GO:0045089      ISS       BP
## 33 ENSG00000154451 GO:0048471      IDA       CC
## 34 ENSG00000154451 GO:0050702      IMP       BP
## 35 ENSG00000154451 GO:0051289      IDA       BP
## 36 ENSG00000154451 GO:0071346      IEP       BP
## 37 ENSG00000154451 GO:0072616      ISS       BP
## 38 ENSG00000154451 GO:1900017      ISS       BP
## 39 ENSG00000154451 GO:1900227      IDA       BP
## 40 ENSG00000264201       <NA>     <NA>     <NA>
## 41 ENSG00000259098       <NA>     <NA>     <NA>
## 42 ENSG00000144369 GO:0016021      IEA       CC
## 43 ENSG00000180658 GO:0003674       ND       MF
## 44 ENSG00000180658 GO:0004930      IEA       MF
## 45 ENSG00000180658 GO:0004984      IEA       MF
## 46 ENSG00000180658 GO:0005886      TAS       CC
## 47 ENSG00000180658 GO:0016021      IEA       CC
## 48 ENSG00000180658 GO:0032154      IDA       CC
## 49 ENSG00000180658 GO:0032467      IMP       BP
## 50 ENSG00000180658 GO:0032956      IMP       BP
## 51 ENSG00000180658 GO:0050911      IEA       BP
## 52 ENSG00000180658 GO:0055037      IDA       CC
## 53 ENSG00000180658 GO:0090543      IDA       CC
## 54 ENSG00000180658 GO:0097431      IDA       CC
## 55 ENSG00000180658 GO:1990023      IDA       CC

goid <- select(org, ensid, "GO", "ENSEMBL") %>% as_tibble()

## 'select()' returned 1:many mapping between keys and columns

head(goid)

## # A tibble: 6 x 4
##   ENSEMBL         GO         EVIDENCE ONTOLOGY
##   <chr>           <chr>      <chr>    <chr>   
## 1 ENSG00000267086 <NA>       <NA>     <NA>    
## 2 ENSG00000184260 GO:0000786 IEA      CC      
## 3 ENSG00000184260 GO:0003674 ND       MF      
## 4 ENSG00000184260 GO:0005634 HDA      CC      
## 5 ENSG00000184260 GO:0005634 IDA      CC      
## 6 ENSG00000184260 GO:0008150 ND       BP

library(GO.db)

## 

columns(GO.db)

## [1] "DEFINITION" "GOID"       "ONTOLOGY"   "TERM"

keytypes(GO.db)

## [1] "DEFINITION" "GOID"       "ONTOLOGY"   "TERM"

gid <- unique(head(goid$GO))
gid

## [1] NA           "GO:0000786" "GO:0003674" "GO:0005634" "GO:0008150"

AnnotationDbi::select(GO.db, gid, "TERM", "GOID")

## 'select()' returned 1:1 mapping between keys and columns

##         GOID               TERM
## 1       <NA>               <NA>
## 2 GO:0000786         nucleosome
## 3 GO:0003674 molecular_function
## 4 GO:0005634            nucleus
## 5 GO:0008150 biological_process

AnnotationHub Cloud resource

Good alternative: annotation hub (web service) –> database

library(AnnotationHub)

beware if there are multiple versions of the same thing (check metadata columns)

hub <- AnnotationHub()

## snapshotDate(): 2019-05-02

query(hub, "^org\\.")

## AnnotationHub with 1710 records
## # snapshotDate(): 2019-05-02 
## # $dataprovider: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
## # $species: Escherichia coli, 'Chlorella vulgaris'_C-169, 'Klebsiella a...
## # $rdataclass: OrgDb
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass,
## #   tags, rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["AH70563"]]' 
## 
##             title                                                     
##   AH70563 | org.Ag.eg.db.sqlite                                       
##   AH70564 | org.At.tair.db.sqlite                                     
##   AH70565 | org.Bt.eg.db.sqlite                                       
##   AH70566 | org.Cf.eg.db.sqlite                                       
##   AH70567 | org.Gg.eg.db.sqlite                                       
##   ...       ...                                                       
##   AH73812 | org.Plasmodium_vivax.eg.sqlite                            
##   AH73813 | org.Burkholderia_mallei_ATCC_23344.eg.sqlite              
##   AH73814 | org.Bacillus_cereus_(strain_ATCC_14579_|_DSM_31).eg.sqlite
##   AH73815 | org.Bacillus_cereus_ATCC_14579.eg.sqlite                  
##   AH73816 | org.Schizosaccharomyces_cryophilus_OY26.eg.sqlite

hub["AH70564"]

## AnnotationHub with 1 record
## # snapshotDate(): 2019-05-02 
## # names(): AH70564
## # $dataprovider: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/
## # $species: Arabidopsis thaliana
## # $rdataclass: OrgDb
## # $rdatadateadded: 2019-04-29
## # $title: org.At.tair.db.sqlite
## # $description: NCBI gene ID based annotations about Arabidopsis thaliana
## # $taxonomyid: 3702
## # $genome: NCBI genomes
## # $sourcetype: NCBI/ensembl
## # $sourceurl: ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/, ftp://ftp.ensembl....
## # $sourcesize: NA
## # $tags: c("NCBI", "Gene", "Annotation") 
## # retrieve record with 'object[["AH70564"]]'

atorg <- hub[["AH70564"]]

## downloading 0 resources

## loading from cache 
##     'AH70564 : 77310'

zoo <- hub["AH70564"]

TxDb packages

Look for available packages

BiocManager::available("^TxDb\\.Hsap") %>% tibble::enframe(name=NULL)

## # A tibble: 5 x 1
##   value                                      
##   <chr>                                      
## 1 TxDb.Hsapiens.BioMart.igis                 
## 2 TxDb.Hsapiens.UCSC.hg18.knownGene          
## 3 TxDb.Hsapiens.UCSC.hg19.knownGene          
## 4 TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts
## 5 TxDb.Hsapiens.UCSC.hg38.knownGene

query(hub, "^TxDb\\.Hsap")

## AnnotationHub with 6 records
## # snapshotDate(): 2019-05-02 
## # $dataprovider: UCSC
## # $species: Homo sapiens
## # $rdataclass: TxDb
## # additional mcols(): taxonomyid, genome, description,
## #   coordinate_1_based, maintainer, rdatadateadded, preparerclass,
## #   tags, rdatapath, sourceurl, sourcetype 
## # retrieve records with, e.g., 'object[["AH52256"]]' 
## 
##             title                                             
##   AH52256 | TxDb.Hsapiens.BioMart.igis.sqlite                 
##   AH52257 | TxDb.Hsapiens.UCSC.hg18.knownGene.sqlite          
##   AH52258 | TxDb.Hsapiens.UCSC.hg19.knownGene.sqlite          
##   AH52259 | TxDb.Hsapiens.UCSC.hg19.lincRNAsTranscripts.sqlite
##   AH52260 | TxDb.Hsapiens.UCSC.hg38.knownGene.sqlite          
##   AH70591 | TxDb.Hsapiens.UCSC.hg38.knownGene.sqlite

library(TxDb.Hsapiens.UCSC.hg38.knownGene)

## Loading required package: GenomicFeatures

txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene

# Look for exons
ex <- exons(txdb)
ex

## GRanges object with 647025 ranges and 1 metadata column:
##                    seqnames        ranges strand |   exon_id
##                       <Rle>     <IRanges>  <Rle> | <integer>
##        [1]             chr1   11869-12227      + |         1
##        [2]             chr1   12010-12057      + |         2
##        [3]             chr1   12179-12227      + |         3
##        [4]             chr1   12613-12697      + |         4
##        [5]             chr1   12613-12721      + |         5
##        ...              ...           ...    ... .       ...
##   [647021] chrUn_GL000220v1 155997-156149      + |    647021
##   [647022] chrUn_KI270442v1 380608-380726      + |    647022
##   [647023] chrUn_KI270442v1 217250-217401      - |    647023
##   [647024] chrUn_KI270744v1   51009-51114      - |    647024
##   [647025] chrUn_KI270750v1 148668-148843      + |    647025
##   -------
##   seqinfo: 595 sequences (1 circular) from hg38 genome

transcripts(txdb)

## GRanges object with 226811 ranges and 2 metadata columns:
##                    seqnames        ranges strand |     tx_id
##                       <Rle>     <IRanges>  <Rle> | <integer>
##        [1]             chr1   11869-14409      + |         1
##        [2]             chr1   12010-13670      + |         2
##        [3]             chr1   29554-31097      + |         3
##        [4]             chr1   30267-31109      + |         4
##        [5]             chr1   30366-30503      + |         5
##        ...              ...           ...    ... .       ...
##   [226807] chrUn_GL000220v1 155997-156149      + |    226807
##   [226808] chrUn_KI270442v1 380608-380726      + |    226808
##   [226809] chrUn_KI270442v1 217250-217401      - |    226809
##   [226810] chrUn_KI270744v1   51009-51114      - |    226810
##   [226811] chrUn_KI270750v1 148668-148843      + |    226811
##                      tx_name
##                  <character>
##        [1] ENST00000456328.2
##        [2] ENST00000450305.2
##        [3] ENST00000473358.1
##        [4] ENST00000469289.1
##        [5] ENST00000607096.1
##        ...               ...
##   [226807] ENST00000619779.1
##   [226808] ENST00000620265.1
##   [226809] ENST00000611690.1
##   [226810] ENST00000616830.1
##   [226811] ENST00000612925.1
##   -------
##   seqinfo: 595 sequences (1 circular) from hg38 genome

exByTx = exonsBy(txdb, "tx")
exByTx

## GRangesList object of length 226811:
## $1 
## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames      ranges strand |   exon_id   exon_name exon_rank
##          <Rle>   <IRanges>  <Rle> | <integer> <character> <integer>
##   [1]     chr1 11869-12227      + |         1        <NA>         1
##   [2]     chr1 12613-12721      + |         5        <NA>         2
##   [3]     chr1 13221-14409      + |         8        <NA>         3
## 
## $2 
## GRanges object with 6 ranges and 3 metadata columns:
##       seqnames      ranges strand | exon_id exon_name exon_rank
##   [1]     chr1 12010-12057      + |       2      <NA>         1
##   [2]     chr1 12179-12227      + |       3      <NA>         2
##   [3]     chr1 12613-12697      + |       4      <NA>         3
##   [4]     chr1 12975-13052      + |       6      <NA>         4
##   [5]     chr1 13221-13374      + |       7      <NA>         5
##   [6]     chr1 13453-13670      + |       9      <NA>         6
## 
## $3 
## GRanges object with 3 ranges and 3 metadata columns:
##       seqnames      ranges strand | exon_id exon_name exon_rank
##   [1]     chr1 29554-30039      + |      10      <NA>         1
##   [2]     chr1 30564-30667      + |      13      <NA>         2
##   [3]     chr1 30976-31097      + |      14      <NA>         3
## 
## ...
## <226808 more elements>
## -------
## seqinfo: 595 sequences (1 circular) from hg38 genome

which.max(lengths(exByTx))

## 31375 
## 31375

Exons from the gene:

exByTx[31375]

## GRangesList object of length 1:
## $31375 
## GRanges object with 363 ranges and 3 metadata columns:
##         seqnames              ranges strand |   exon_id   exon_name
##            <Rle>           <IRanges>  <Rle> | <integer> <character>
##     [1]     chr2 178807212-178807423      - |     92853        <NA>
##     [2]     chr2 178804552-178804655      - |     92850        <NA>
##     [3]     chr2 178802138-178802341      - |     92847        <NA>
##     [4]     chr2 178800395-178800682      - |     92846        <NA>
##     [5]     chr2 178799825-178799910      - |     92845        <NA>
##     ...      ...                 ...    ... .       ...         ...
##   [359]     chr2 178529960-178530116      - |     92476        <NA>
##   [360]     chr2 178528528-178529219      - |     92475        <NA>
##   [361]     chr2 178528274-178528427      - |     92474        <NA>
##   [362]     chr2 178527446-178527748      - |     92473        <NA>
##   [363]     chr2 178525989-178527307      - |     92471        <NA>
##         exon_rank
##         <integer>
##     [1]         1
##     [2]         2
##     [3]         3
##     [4]         4
##     [5]         5
##     ...       ...
##   [359]       359
##   [360]       360
##   [361]       361
##   [362]       362
##   [363]       363
## 
## -------
## seqinfo: 595 sequences (1 circular) from hg38 genome

TCGA

BiocManager::available("curate")

## [1] "curatedAdipoChIP"       "curatedAdipoRNA"       
## [3] "curatedBladderData"     "curatedBreastData"     
## [5] "curatedCRCData"         "curatedMetagenomicData"
## [7] "curatedOvarianData"     "curatedTCGAData"

Download data from TCGA

library(curatedTCGAData)

## Loading required package: MultiAssayExperiment

curatedTCGAData()

## Please see the list below for available cohorts and assays

## Available Cancer codes:
##  ACC BLCA BRCA CESC CHOL COAD DLBC ESCA GBM HNSC KICH
##  KIRC KIRP LAML LGG LIHC LUAD LUSC MESO OV PAAD PCPG
##  PRAD READ SARC SKCM STAD TGCT THCA THYM UCEC UCS UVM 
## Available Data Types:
##  CNACGH CNACGH_CGH_hg_244a
##  CNACGH_CGH_hg_415k_g4124a CNASeq CNASNP
##  CNVSNP GISTIC_AllByGene GISTIC_Peaks
##  GISTIC_ThresholdedByGene Methylation
##  Methylation_methyl27 Methylation_methyl450
##  miRNAArray miRNASeqGene mRNAArray
##  mRNAArray_huex mRNAArray_TX_g4502a
##  mRNAArray_TX_g4502a_1
##  mRNAArray_TX_ht_hg_u133a Mutation
##  RNASeq2GeneNorm RNASeqGene RPPAArray

curatedTCGAData("BRCA")

##                                         Title DispatchClass
## 31                       BRCA_CNASeq-20160128           Rda
## 32                       BRCA_CNASNP-20160128           Rda
## 33                       BRCA_CNVSNP-20160128           Rda
## 35             BRCA_GISTIC_AllByGene-20160128           Rda
## 36                 BRCA_GISTIC_Peaks-20160128           Rda
## 37     BRCA_GISTIC_ThresholdedByGene-20160128           Rda
## 39  BRCA_Methylation_methyl27-20160128_assays        H5File
## 40      BRCA_Methylation_methyl27-20160128_se           Rds
## 41 BRCA_Methylation_methyl450-20160128_assays        H5File
## 42     BRCA_Methylation_methyl450-20160128_se           Rds
## 43                 BRCA_miRNASeqGene-20160128           Rda
## 44                    BRCA_mRNAArray-20160128           Rda
## 45                     BRCA_Mutation-20160128           Rda
## 46              BRCA_RNASeq2GeneNorm-20160128           Rda
## 47                   BRCA_RNASeqGene-20160128           Rda
## 48                    BRCA_RPPAArray-20160128           Rda

mae <- curatedTCGAData("BRCA", c("GISTIC_AllByGene", "RNASeq2GeneNorm"), dry.run=FALSE)

## snapshotDate(): 2019-04-29

## see ?curatedTCGAData and browseVignettes('curatedTCGAData') for documentation

## downloading 0 resources

## loading from cache 
##     'EH588 : 588'

## see ?curatedTCGAData and browseVignettes('curatedTCGAData') for documentation

## downloading 0 resources

## loading from cache 
##     'EH596 : 596'

## see ?curatedTCGAData and browseVignettes('curatedTCGAData') for documentation

## downloading 0 resources

## loading from cache 
##     'EH587 : 587'

## see ?curatedTCGAData and browseVignettes('curatedTCGAData') for documentation

## downloading 0 resources

## loading from cache 
##     'EH590 : 590'

## see ?curatedTCGAData and browseVignettes('curatedTCGAData') for documentation

## downloading 0 resources

## loading from cache 
##     'EH599 : 599'

## harmonizing input:
##   removing 12081 sampleMap rows not in names(experiments)

mae

## A MultiAssayExperiment object of 2 listed
##  experiments with user-defined names and respective classes. 
##  Containing an ExperimentList class object of length 2: 
##  [1] BRCA_GISTIC_AllByGene-20160128: SummarizedExperiment with 24776 rows and 1080 columns 
##  [2] BRCA_RNASeq2GeneNorm-20160128: SummarizedExperiment with 20501 rows and 1212 columns 
## Features: 
##  experiments() - obtain the ExperimentList instance 
##  colData() - the primary/phenotype DataFrame 
##  sampleMap() - the sample availability DataFrame 
##  `$`, `[`, `[[` - extract colData columns, subset, or experiment 
##  *Format() - convert into a long or wide DataFrame 
##  assays() - convert ExperimentList to a SimpleList of matrices

names(assays(mae))

## [1] "BRCA_GISTIC_AllByGene-20160128" "BRCA_RNASeq2GeneNorm-20160128"

mae[["BRCA_RNASeq2GeneNorm-20160128"]]

## class: SummarizedExperiment 
## dim: 20501 1212 
## metadata(0):
## assays(1): ''
## rownames(20501): A1BG A1CF ... psiTPTE22 tAKR
## rowData names(0):
## colnames(1212): TCGA-3C-AAAU-01A-11R-A41B-07
##   TCGA-3C-AALI-01A-11R-A41B-07 ... TCGA-Z7-A8R5-01A-42R-A41B-07
##   TCGA-Z7-A8R6-01A-11R-A41B-07
## colData names(0):

example: calculate distribution of library sizes

mae[["BRCA_RNASeq2GeneNorm-20160128"]] %>% 
  assay () %>%
  colSums() %>% 
  hist()